Pattern formation of granular avalanches with vortex convection 
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It has been known that the granular flow of polystyrene particles down a slope forms a wavy 
pattern with many heads at the moving front of the resulting avalanche. In experiments of gran- 
ular flow using low-density particles, the instability of the moving front and the subsequent head 
formation are driven by gravity and air drag. To elucidate the instability mechanism of granular 
avalanches, we propose a particle method considering gravity as the driving force for the avalanche, 
the contact interaction between granular particles, and the long-range interaction between granular 
particles through the ambient fluid as a type of air drag. Using this model, we simulate the head 
formation at the moving front of the avalanche, and we investigate the particle flow caused by the 
air drag. It is found that the air drag destabilizes the shape of the avalanche that deforms into a 
wavy pattern, leading to the generation of a pair of granular vortex convection currents inside the 
head. Further, the relationship between the particle radius and head size is found to satisfy the 
positive linear scaling law. Moreover, the hydrodynamic interaction between the particles causes 
their aggregation at the moving front of the avalanche, and this aggregation effect generates the 
head-tail structure. These numerical results are qualitatively consistent with the results of previous 
experiments. 
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Avalanches, generally taken as a class of massive land- 
slide phenomena, involve gravity currents and density 
currents, for instance, snow avalanches, debris flows, and 
pyroclastic flows. These flows slide down a slope as a 
mixture of solid and fluid and form several kinds of typ- 
ical structures, one of which is termed as the head-tail 
structure Concretely, at the moving front of an 

avalanche, a large head is formed through the clustering 
of large particles, whereas at the rear end, an elongated 
tail is formed by the trailing smaller particles. Among 
potential factors for the formation of the head-tail struc- 
ture, the ratio of air drag to gravity is considered as 
the most relevant. For example, a granular flow consist- 
ing of low-density particles such as polystyrene particles 
forms the head-tail structure as well as the wavy pattern 
with many heads at the moving front of the avalanche 
In addition, these experiments show that the pos- 
itive linear relation between head size and particle ra- 
dius is a rough approximation. In contrast, experiments 
using high-density particles such as glass beads do not 
generate the head-tail structure, although they form a 
wavy pattern similar to that obtained when using low- 
density particles @- This difference in the formation of 
the avalanche structure due to the difference in particle 
density can be estimated based on whether or not the air 
drag efficiently acts on the particles [9] . 

In order to explain these facts by theoretical means, 
several models have recently been proposed; the fluid 
flow model assumes an avalanching body as a mass of 
fluid, the mass center model assumes an avalanche as a 
huge single particle, and so on [Iol - |T3 |. However, the 
materials constituting avalanches that are used in these 
models are granular materials such as polystyrene or glass 
beads, which are definitely not fluid. Moreover, the sin- 
gle particle model does not treat the interaction between 
granular particles, which is speculated to play a crucial 
role in the dynamics of an avalanche. To overcome the 



above shortcomings these models, we employ a particle 
method considering the contact between granular parti- 
cles and the fluid effect as the air drag. In this study, we 
investigate the head formation in a granular avalanche 
by considering the particle radius, and we clarify the in- 
stability mechanism of granular avalanches, which leads 
to the generation of the head. 

In our model, a granular flow slides down a slope that 
is steeper than the angle of repose, and this model is 
roughly based on three basic assumptions: First, the 
granular flow consisting of spherical (three-dimensional) 
particles only moves along either a 2D inclined slope or a 
2D vertical cross-section extending over a slope. Second, 
only the translational motion of particles along either of 
the two 2D spaces is considered, that is, the rotational 
motion of particles is ignored. Third, as regards the main 
force acting on the particles, we consider three types of 
forces: (i) gravity as the dominant driving force for the 
avalanche, (ii) the repulsive force between particles that 
is caused by the excluded volume effect, and (iii) the drag 
force due to the fluid. The specific form of each force is 
described as below. 

(i) Gravity: Considering the density difference between 
the granular particles and the fluid, the gravity of the zth 
particle F g i is expressed as 

F gi = -e v Vig(p p - pf), 

where e y denotes the unit vector parallel to the y-axis, Vi 
denotes the volume of the ith particle, g denotes the grav- 
itational acceleration, and p v and pf denotes the densities 
of the particle and the fluid, respectively. 

(ii) Repulsive force: We ignore the rotational motion 
of the particles, and we assume that the ith particle ex- 
periences only a normal directional force F r j, which acts 
as the repulsive force caused by contact with the jth par- 
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ticle and is represented as the elastic spring force 
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where denotes the relative angle of r = — to 
the positive direction of the x-axis, where r 4 and Vj rep- 
resent the respective central coordinates of the ith and 
jth particles, k n denotes the spring constant, and 5 de- 
notes the overlap distance between two particles. Here, 
5 is expressed in terms of the particle radius a and the 
interparticle distance r — [xi — Xj) 2 + (jji — yj) 2 as 



+ a,j — r (r < a.- L + a,j ; contact) 
(else; noncontact) . 



(iii) Drag force by fluid: The long-range interaction 
between the particles though the fluid is experienced by 
individual particles as a type of air drag force. Con- 
sequently, we assume that the motion of each particle 
immediately propagates though the fluid, and we use 
the Rotne-Prager tensor J(r) [l5| to denote the effec- 
tive long-range interaction between particles. This ten- 
sor is the exact solution of the Stokes equation in cases 
in which the fluid contains spherical particles, and the 
specific form is expressed as 



J(r) = - 

r 
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where I denotes the unit tensor. Here, Eq. ([!} indicates 
the fluid effect caused by the jth particle on the ith par- 
ticle, and rr in Eq. (JlJ denotes the matrix in which each 
element corresponds to the tensor J; for instance, the ten- 
sor element 3 xy is expressed as (rr)^ = (xi—Xj)(yi — yj). 

Moreover, by ignoring the particle inertia and only con- 
sidering the three abovementioned effects, the ith particle 
velocity v,; is determined as 
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where F* = F g i + F r ,, p denotes the viscosity coefficient 
of the fluid, and Uj(j) denotes the induced velocity gen- 
erated by the jth particle. Further, the first term in Eq. 
@ is the expression for Stokes' drag formula. 

Numerical simulations of the particle model are con- 
ducted on a flat slope with a constant inclination an- 
gle (9 — 45°), where particles are confined in two differ- 
ent kinds of two-dimensional planes depending on model 
(Fig. Q}: model-I) particles move along the 2D surface 
of the slope, and model-II) particles move within the 2D 
vertical cross-section extending over the slope. It is to be 
noted that model-I addresses the instability to generate 
wavy pattern at the moving front of an avalanche from 
the top-view perspective, where the effect of the inclina- 
tion angle is reflected in the gravity term as g sin 9 and 
the friction between the particle and the slope is ignored. 



In contrast, model-II focuses on the head-tail structure 
formation and assumes that the collision of particles with 
the slope is described as the elastic spring force in a man- 
ner similar to the force involved in the collision between 
particles. 

As regards the specific materials of the granular par- 
ticle and fluid, we select a low-density particle such as 
polystyrene with the choice of the fluid as air for the pur- 
pose of experimental comparison. Concretely, the density 
p p and the spring constant k n of the granular particle are 
given as p p = 20 kg m -3 and k n — 10 N m _1 , respec- 
tively, whereas the density pf and the viscosity coeffi- 
cient p of the fluid are fixed as pf = 1.2 kg m~ 3 and 
p = 1.82 x 10~ 5 Pa s, respectively. As the control pa- 
rameter, we vary only the particle radius a p in the range 
of the order of millimeters, that is, 

a p = (1 ± 0.05)s x 10~ 3 m 
s = {1.0, 2.5, 5.0} (small, medium, large), (3) 

where s denotes the criteria radius. Each particle used 
for the calculations is randomly selected within the range 
[0.95 x s, 1.05 x s] to avoid crystallization of the particles. 
Further, both simulations of model-I and II are carried 
out 20 times with N = 2000 particles for different sets of 
particle radius and initial positions. 

In model-I, we use two different configurations of the 
initial particle setup: i) circular setup (Fig. HJa)) and 
ii) rectangular setup (Fig. G2a)), which are introduced 
to focus on the formation of a single head and the struc- 
tural instability at the moving front of avalanche, respec- 
tively. The spatial dimensions of both initial setups are 
chosen to be sufficiently large to pack N = 2000 par- 
ticles into each setup depending on the criteria radius. 
Hence, the radius of the circular initial setup R and the 
lateral length of the rectangular initial setup L are set as 
R = {0.065,0.17,0.34} and L = {0.35,0.85,1.7} (small, 
medium, large), respectively. For the directional length 
of the slope of the rectangular initial setup, we fix the 
initial length as 10 particles. It is to be noted that the 
particles are randomly placed inside each setup, avoid- 
ing overlap between particles, and the packing fraction is 
nearly 0.5. In addition, the particle velocity v is derived 
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FIG. 1. Two different kinds of two-dimensional planes: (a) 
Model-I (top view), (b) Model-II (side cross-sectional view). 
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FIG. 2. Head formation for circular initial setup with 
s = 5.0; (a) t = s. (b) Particle velocity field against migra- 
tion speed of avalanche v a at t — 2.45 x 10 -3 . The y com- 
ponent of the velocity vector is described as v y — v a , where 
v y denotes the y component of the particle velocity, (c) Def- 
inition of width and layer thickness in head. The gray and 
black circles indicate the constituents of head and targets se- 
lected to measure head size, respectively, (d) Time evolutions 
of head size rescaled by criteria diameter. 



from the particle model with a single particle as 
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(4) 



Equation (0| indicates that the particle radius a p strongly 
governs the particle velocity, and thus, numerical sim- 
ulations are run until the simulation time reaches t = 
2 x 10~ 2 , 10~ 2 , 5 x 10~ 3 s (small, medium, large) depend- 
ing on the criteria radius (Eq. (J3|)). 

Next, we begin with the study of the head forma- 
tion process using model-I with a circular initial setup. 
Initially, the particles are uniformly located within this 
setup (Fig. [U[a)), and they aggregate around the mov- 
ing front as they travel down the slope, as described by 
Eq. (0); eventually, the aggregated particles form a sin- 
gle head. To observe the relative motion of particles at 
the head, we visualize the particles velocity relative to 
the mean migration speed of the head during the head 
formation. Figure [U[b) shows that the particles along 
the rear (upper) verge of the densely aggregated band 
of the head are pulled forward by the neighboring parti- 
cles and move toward the center of the head, whereas the 
particles along the front (lower) verge of the head move 
away from the center of the head. In this manner, the 
granular avalanche migrates downward while a pair of 
granular vortices is maintained inside the head, through 
which process the head shrinks with time by leaving be- 
hind a certain portion of particles. It is to be noted that 



similar processes are observed independent of the criteria 
radius (Eq. ©). 

In order to measure the head size, we define the width 
w and layer thickness a of the head as follows: First, to 
discern the head from the other structural parts of an 
avalanche, we choose particles according to the quantity 
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where S r = 4 x s denotes the search range. Here, nl i in- 
dicates that the ith particle includes indirect neighbors 
through the jth particle, that is, n' yields the quasilo- 
cal particle number density. If the zth particle satisfies 
the condition v! i > n' max /4, we select this particle as a 
constituent of the head, which is indicated as a gray (or 
black) circle in Fig. HJc). Furthermore, we choose four 
target (i.e., black) particles within the head, which are 
located at the foremost position in the head, the fore- 
front along the rear verge of the head, and the outmost 
positions in right and left arms of the head. According 
to the locations of these four black circles, we define the 
size of the head; the head width w is defined as the lat- 
eral length between the two outmost particles, whereas 
the head layer thickness 8 is the slope directional length 
from the foremost particle to that along the rear verge in 
the head. The plot in Fig. H^d) indicates that the head 
width shrinks with time to reach a finite width of about 
20 particles during the simulation time independent of 
the criteria radius (Eq. ([3])). Moreover, we found a 
strong positive correlation (correlation coefficient: 0.956) 
between the width w = w/(2s) and the layer thickness 
8 = 8/ '(2s) rescaled by the criteria diameter. 

Next, the simulations using model-I with a rectangular 
initial setup show a time evolution in which the air drag 
destabilizes the initially straight front of the avalanche 
to deform into a wavy pattern with many heads (Figs. 
|3ja) and (b)). For this growth process of instability, 
we speculate that the dominant factor is the long-range 
fluid-mediated interaction between particles according to 
Eq. ([1]). This is because; the section of particles located 
around the middle height in the rectangular initial setup 
consist of a larger number of fluid-mediated interaction 
particles than those located near the verges of the setup; 
thus, the velocities of the former particles are larger than 
those of the latter, as reflected by the second term in Eq. 
@. This velocity difference between the middle, frontal 
and rear particles gives rise to the aggregation of parti- 
cles at the moving front of the avalanche as in the case 
of the circular initial setup. However, in this case of the 
laterally long initial setup, the randomness of the par- 
ticle radius and initial inhomogeneity contributes to the 
subsequent evolution of the lateral inhomogeneity of the 
frontal pattern of the avalanche, which breaks down into 
a multi-head structure (Fig. GUb)). 

The heads in the wavy structure are characterized 
by a high particle number density, and we determine 
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FIG. 3. Instability at moving front of avalanche for rectangu- 
lar initial setup with s — 5.0; (a) t = s and (b) t = 3x 10 -3 
s. (c) Scatter diagram for particle number of head and head 
width. The solid line per criteria radius is obtained by the 
least square method, (d) Average width against particle ra- 
dius. The vertical bars indicate absolute deviations, whereas 
the line passing through the origin is also fitted by the least 
square. 



each head size according to the previously introduced ap- 
proach (Fig. Ell a))- It is to be noted that numerical simu- 
lations with the circular initial setup were conducted for 
different time intervals depending on the criteria radius 
(Eq. ([3])) such that the head sizes rescaled by the crite- 
ria diameter reached a certain value ( w ' ~ 20) during the 
simulation time (Fig. [2jd)) . Further, in the second case 
of the simulations with the rectangular initial setup, the 
time interval of the simulations is set depending on the 
criteria radius such that the head formation is clearly ob- 
served, where the head size is measured at the end of the 
calculation. The mean numbers of emergent heads per 
simulation are 17.2, 13.95, and 14.85 (for small, medium, 
and large particles). The results in Fig. [3jc) indicate 
that the width of the head increases with the number 
of particles constituting the head (i.e. gray or black cir- 
cles in Fig. [2jc) ) , and the relation between these two 
parameters can be roughly fitted by straight lines. To 
study the relation between the particle radius and the 
head size, we calculate the average width and obtain a 
linear relationship between the particle radius and the 
head width. These results quantitatively correspond to 
a previous experiment performed using polystyrene, and 
they show a similarity to the granular Rayleigh-Taylor 
instability observed in sedimenting suspensions fl6l . 
Finally, we establish whether model-II describing the 



granular avalanche in a vertical cross-section can re- 
produce the intrinsic head-tail structure for low-density 
granular particles. As the initial setup, the packing area 
of particles is prepared by a stopper extending perpen- 
dicular to the slope with inclination angle 6 — 45° (Fig. 
Eta)). The size of this packing area is given according to 
the vertical depth D, which is set as D = {0.115, 0.3, 0.6} 
(small, medium, large) where the packing fraction is 
nearly 0.5. The numerical simulations show that the 
particles converge to the moving front of the avalanche 
via the long-range fluid-mediated interaction according 
to Eq. ((T|), and subsequently, the thickness of the front 
part becomes more than that of the rear part, that is, 
the granular avalanche forms the head-tail structure (Fig. 
BJb)). Moreover, the head shrinks with time while retain- 
ing its structure because the particles located at the top 
of the head move backward. 

In this study, we employed a particle method consider- 
ing the fluid-mediated non-local particle-particle interac- 
tion using the Rotne-Prager tensor, and we investigated 
the pattern formation of granular avalanches in two dif- 
ferent types of 2D space: (i) surface of constant slope and 
(ii) vertical cross-section extending over the slope. The 
numerical simulations using (i) showed that the fluid- 
mediated particle-particle interaction deforms the gran- 
ular avalanches into a wavy pattern with many heads at 
the moving front. In addition, each head migrated down- 
ward while maintaining a granular vortex convection, for 
which the relation among particle radius, head thickness, 
and head width satisfied the linear scaling law. Moreover, 
the simulations using (ii) were successful in the reproduc- 
tion of a characteristic head-tail structure independent of 
particle radius. These results qualitatively correspond to 
the results of previous experiments. 




FIG. 4. Head-tail structure of avalanche in model-II with 
s = 5.0. (a) t = s. D denotes the vertical depth of the 
packing area, (b) Particle velocity along slope at t = 2 x 10~ 3 
s. The particles with velocities between the migration speed 
of avalanche v a and maximum velocity u ma x are characterized 
by the black color gradation. 
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